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ABSTRACT 

A method is presented for systematic and geodetic correction data calculation. 

It is based on presentation of image distortions as a sum of nominal distortions 
and linear effects caused by variation of the spacecraft position and attitude 
variables from their nominals. The method may be used for both MSS and TM image 
data and it is incorporated into the processing by means of mostly offline 
calculations. Modeling shows that the maximal errors of the method are of the 
order of 5m at the worst point in a frame; the standard deviations of the average 
errors less than .8m. 

INTRODUCTION 

The geometric correction of the Landsat-type imagery typically proceeds in 
two steps. The first, called the Systematic Correction, removes internal distortions 
imported in the raw image data by the sensor mechanism, spacecraft motion, inaccurate 
sensor pointing, earth’s rotation, etc. These partly corrected images still contain 
distortions due to uncertainties in spacecraft position and orientation. The second 
step. Geodetic Correction, removes these residual distortions using refined values 
of the attitude and ephemeris estimates. The refined attitude and ephemeris are 
obtained by filtering of image dislocations at Control Points. 

Application of the geometric correction requires the generation of the Correction 
Data - Systematic (SCD) or Geodetic (GCD) , depending upon the processing step. 

This data is developed on a rectangular grid in input (pixel, scan line) coordinates 
and express the relationship between the input and output map coordinates, within 
a standard World Reference System (WRS) frame. 

The central part of the SCD/GCD generation is the computation of the coordinates of 
the intersection of the sensor’s llne-of-slght vector, with the Earth’s surface 
(lookpoint coordinates). The lookpoint coordinates must then be converted to 
geodetic coordinates followed by mapping into user’s map projection. There are two 
user’s map projections: Space Oblique Mercator (SOM) and either Universe Transverse 

Mercator (UTM) or Polar Stereographic (PS) . 


* Work performed under National Aeronautics and Space Administration 
Contract No. NAS 5-25300. 
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Finally, the data, computed for integer values of pixels and lines, is inter- 
polated to integer values of output map coordinates. The grid spacing is chosen 
so that the data, together with properly defined interpolation techniques, represent 
the output coordinates to the desired precision everywhere in the frame. 

It should be noted that all the calculations are performed twice at each grid 
point, once for each SCD and GCD. They consume a significant amount of the processing 
time, which needs to be minimized. At the same time, there are no essential 
differences between SCD and GCD. Both establish a pointwise transformation, which 
may be written generically as 

^m ^ F(pixel,line,p) , 

where coordinates of a grid point and p is a vector of 

variables characterizing the spacecraft motion, attitude pointing, sensor’s 
mechanism, etc. 


Letting p = p^ + 6 , the sum of nominal values of the variables and the deviation 


from the nominals, in the first approximation 

' n 


X = X + y 5 , 
m m 


( 1 ) 


where X^ are the nominal map coordinates and y is the matrix of the partial 
derivatives (PD) 


y = 


[■f] 


Thus, SCD and GCD may be represented as a sum of the nominal correction data and 
pointwise adjustments. 

This has significant advantages: 

1) It provides a uniform approach to the SCD and GCD computations, considering 
each as one transformation, and 

2) Because the nominal spacecraft motion is known for every WRS frame, the nominal 
coordinates and the partial derivatives may be computed and stored in a Data Base. 


The implementation of such an approach depends a great deal on both the choice 
of an output map projection and 6. An analytic form of mapping not only has to 
allow derivation of the coefficients y, but it should also afford rapid and precise 
online inversion to geodetic coordinates, from which the final map projection can 
be generated. In addition, it is desirable to have the nominal coordinates and the 
partial derivatives, as far as possible, insensitive to global position of the 
frame. Thus, although out technique may be applied to most standard map projections 
(such as UTM or PS), a special intermediate projection, Local Space Oblique Mercator 
(LSOM) , has been employed. The LSOM is the Mercator projection for the sphere, with^ 
local ’equator’ along the nominal spacecraft inertial velocity at the frame center. 
In that projection X and y are longitude-undependent and thus, can be stored only 
for one path. A natural choice of variables 6 is the along-track, cross-track and 
radial deviations in spacecraft position, together with deviations in the attitude 
angles. The nominal spacecraft motion within a frame is assumed to be in a perfect 
circular orbit passing through the frame center. 
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NOTATIONS 


3p (^pl > ®p2 ’ • • • ’ 


\\h\ 


X 

s 



(X 


ml’ 


X 


m2 


) 


n 


~ vector (Kxl matrix) 

- the Euclidean norm of a 

P 

~ transposed matrics A 

spacecraft position vector in earth- 
centered earth-fixed coordinates 

spacecraft position vector in nominal 
spacecraft coordinates 


- output map coordinates 


pointing vector in body coordinates 


g 


pointing vector in local vertical 
spacecraft coordinates 


f 


pointing vector in earth-centered fixed 
coordinates 



coordinates of a lookpoint on the earth 
surface 


h 


V 


distance from spacecraft to earth 
lookpoint 



normalized lookpoint vector. 


R 


- local earth radius at WRS frame center 




earth rotation rate 



earth equatorial and polar radii 




1 , 


®3 = 


agb 


-1 

e 
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deviation in the pitch 

deviation in the roll 

deviation in the yaw 

along track angular deviation 

cross track angular deviation 

relative departure in the radial direction 

time deviation 

(2x1) matrix of the partial derivatives 

of X with respect to 6^ . 
m ^ k 


li = (PpM 



a 

S 


X 



the ^ equivalent ' pitch and roll 
spacecraft orbital rate 
geodetic longitude 
geocentric latitude 
geodetic latitude 


ROTj^(y) = 




ROT^Cil;) = 


10 0 
0 cosip -sini|) 


0 sinij; 

cosip j 

cosip 0 

sinip 

0 1 

0 

-sinip 0 

cosip 

cosip 

-sinip 0 

sinip 

cosip 0 


0 0 1 
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JL ROT.W =T. 

r. 0 0 0 

^10 0 -1 
0 10 

0 0 1 

T„ = 0 0 0 

^ -1 0 0 

0 -1 0 

T„ = 1 0 0 

■^ 0 0 0 


I - the three dimensional identity matrix 

The upper index n indicates the nominal value of a vector. 

T ^ “ active scan time 
act 

T - mirror turnaround time 

round 


THE NOMINAL SCD 
Coordinate Transformations 

The local (instantaneous) spacecraft coordinates are described in terms of 
the unit vectors where points towards the Earth center, the 

vector is along the orbital angular momentum, and ^2“ To ^ Ti is roughly along 
the velocity direction. The local spacecraft coordinates at ^ the WRS center i^ 
called the nominal snacecraft coordinates. The matrix A transforms a vector X in 
earth-centered inertial coordinates to the vector X in nominal spacecraft ^ 
coordinates: ® 

X = M (2) 

s o ■ ^ 

The inertial to earth-centered earth-fixed coordinate transformation is 
defined as 

X ^ e\, (3) 

where E == E • R0T^(f2t). 

o 3 

The matrix E^ gives the time- independent component of the transformation, 

ROT^(f^t) describes the rotation about the earth axis at the rate ^2. We assume that 
t = 0 at the frame center. The corresponding nominal spacecraft to earth fixed 
coordinate transformation may be written as 

— T T— T — 

X = E A X = P X (4) 

s s ^ 


where P = AE. 


8-5 



The unit vector g”, given in body (sensor) coordinates, is transformed to local 
spacecraft coordinates as 

g = ROT (-e ) • ROT (-0J • ROT(-0_) i” (5) 

J y ; Z 17 p 

where 0^, 0^, 0^ are the yaw, roll, and pitch. 

In the nominal spacecraft coordinates, g may be expressed as 

s 

is = ROTj^(y)i 

where y is the angle in the orbit plane between the spacecraft and the f_rame center. 
In the nominal spacecraft motion cos y = - XgS/ |Xs|| , slnY = Xg£/ || ^sl ’ 
and thus, the matrix ROT^^( Y) is known completely. 


A vector X = (X . , X „) in LSOM coordinates is defined as 
m ml’ m2 


X 


‘ml 


R ^ 1 + sin 3 

2 1 - sin 3 


(7) 


X „ = Ra 

m2 


where R is the earth local radius at the frame center. The local polar angles, 
a and 3, are given by 


= sin3 

W 2 = cos3 • sina (8) 

= -cos3 • cosot 


where _ _~1 _ 

W = II ull Au (9) 

and u = (u^,u ,u ) are earth fixed coordinates of the corresponding point on the 
ground. 
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Generation of the nominal SCD 

The nominal coordinates, X^, are computed on a grid, consisting of 2ni+l 
fictitious scan lines, each line containing 2n2+l points. Because the TM scans in 
two directions, it requires two sets of the nominal coordinates, for forward and 
backward scans. The computations may be fulfilled in the following order. 


1. Generate the time of (i,j) point 


Here 


'ij 


T 

-s. c . en g. 


(i-ni-1) 


^act 

2n2 


(j-n2-l) + AT 


T = (T ^ + T ,) 

scene act round 


(N - K.) 
scan 1 


where T is the active scan time, T a turnaround time, N is the 

actual number of scans, and = 1 for^MSS and 2 for TM. scan 

The parameter T = T for backward scans of the TM and zero otherwise. 

scene 

2, Generate 2n2+l unit line-of-sight vectors in body coordinates. An actual 
mirror velocity profile, together with constant sensor’s misalignments may be employed 


3. Compute the spacecraft position vector X and the matrix P at t . 

s Ij 

4. Compute g^ according to (6) . 

5. Transform and g^ in earth fixed coordinates, obtaining the vectors X and f, 
respectively. 

6. Determine the lookpoint coordinates, u ^ , and h from the equations 

u = X + hf (10) 

(u" + u") + u^ b^2 = 1 (11) 

7. Transform u into LSOM coordinates using (9), (8), and (7). 

It is convenient to have all distances in units of the nominal orbit radius. 


THE PARTIAL DERIVATIVES 
Position and Pointing Vectors 

~'n , 

Let be a nominal spacecraft position vector at time t , 6, and 6 be the 

angular along-track and cross-track deviations in spacecraft posihon, and 5^ be 

a relatlv^e deviation in the radial direc_tlon. Then the actual spacecraft position 

vector, X , may be obtained by rotating X’^ through 6, and 6 . This is followed by 

stretching according to the ratio 1-6,:^ ^ ^ 

6 

Xg = R 0 T 2 ( 63 > • ROT^(d^) • if (1-6^) 


8-7 



Similarly, if and 6^ are deviations in the pitch, roll, and yaw, the actual 

(unit) pointing vector should be written as 

ig= ROT 2 ( 65 ) • R 0 T^( 6 g)i^ , 
where is defined by (5) and (6) as 

i" = ROT^(y) R0T2(-63)R0T2(-52) * R0T^(-5^) 

Let t = t + 6- . Remembering that P = P(t ) = AE * ROT„(fit ) , we can write 

Q / O J 

T 

P at time t as 

p'^(t) = ROT^ (-n(t^+5^))E^V = R0T^(-a6^)' 

•ROT„(-J^t )e'^a'^ = R0T„(-a6 )P'^ 

3 o o i / 


and the actual position and pointing vectors in earth fixed coordinates as 
X = R0T2(-J^67)P^R0T2(§5)R0T^( 5^)X^(l-6g) 
f = R0T^(-!:26^)p\0T2(63)R0Tj^(d^)R0T^(Y) ■. 

R 0 T^(- 63 > • R 0 T 2 (- 62 ) ‘ R0T^(-6^)i’^ 

The linear terms of the Taylor series expansions of X and f in the vicinity of 
6 ^ = 0 (i = 1,2,. . . ,7) give 

X = p’^ (I + 1 ^ 5 ^+ 1265 - 5 ^-^ 136 ^)^ 

r ^ 

f = P^ I^I+T^ROTj^Cy) 6 ^+T2R0T^(y) 55-ROT^ (y) • '^i'^i" 

rot^(y) 5J 

Here we used the fact that 

— ROT, (ijj) , ^ = T. 
dtjj 1 ^ rp=0 1 
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Introducing = p'^ROTj(Y)g^ and 

» we finally have 

X = X" + 


f = + p'^(T^ 6 ^+T 26 ^-OT 36 ^)Pf"- 


( 12 ) 


- P\0T (y)(’5^ T 6 .)g^ 
^ i=l ^ ^ 


The Partial Derivatives of Lookpoint Coordinates 

Henceforth, we will use a prime to denote the matrix of PD with respect to 5 
computed at the nominal point. From (10) it follows that 


M 


9u 

86 


-1 —1 -1 1 - 
= u = X + hf + h f 


(13) 


Introducing e^^ = 62 = 1 and e^ = ^e^e ^ ™ay be rewritten as 


ef (x'.+hf.) = a 

^ XIX 

i=l 


or. 


'( 5“ f.ef) + 2h(T“f,X.ef) + ]^X? e? = a 

/L^ X X XXX ^ — X X 


2 

i 

e 


Differentiating the last expression as a inpllclt function of h gives 

2 


^ (Xi + hfh (X■^ + hf.-)e? 
. _ Zr u,ei(Xj+hf}) 

r\ 

f,u,e 


i^Ci 


1 


and, after substitution h in (13), we have 

r 

l^k 


where 


\_(x} + h4) ^ 


(i,k = 1,2,3) 
Co = -^fiUjei 
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Using the matrix notations, may be expressed as 
= C(X^ + hf^) J 


( 14 ) 


where C is a 3x3 matrix with the elements 
C, 


. = 1 + U, f, e^C - 
kk k k k o 


C, . = u.f 
kj J k j o 


(15) 


Transformation to LSOM Coordinates 

To transform the lookpoint coordinates_^ u,_^to the LSOM coordinates, they 
must be represented in the normalized form V = u* . 


Differentiation of V Yields 


V. 


= (U^ “ \\^\^ 3 ^ U-u!") (lull ^ and introducing the matrix B with 

K. K. K. £ I 1 X 


the elements 


\k - 


(16) 


IW9 


-2 


k j 


may be written as 


- Hull ^ BU^ = Ijull ^ BC(X^ + hf^) 


-1 


(17) 


The next step is transformation of V to W and then to X^^^. From (9) and (17) 


it follows that 


= 


-1 


ABC(X^ + hf^) 


(18) 


From (7) and (8) it follows that 


1 1 + Wi 

•ml' 


( 19 ) 


X ^ = R*arctan (-W^/W„) 
mz z J 


and therefore, 


1 2-1 1 

X , = ,R(1-WJ 
mi i i 


X*2 - R(1-w2)(W2w‘ - W3WI). 
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Introducing the matrix 


D = /I 0 0 

0 -Wj Wj 


( 20 ) 


and using (18), we have 


—1 2 - 1—1 

X = R(i-w:) w = 
m i 


R jju II "^(1-W^) ^ DABC(X^ + hfb 


( 21 ) 


obtain the final result, we must substitude an explicit expression for 
X + hf , which follows immediately from (12) : 


-I 


'“-hP^OT^Cy) g’^ 

k= 1,2,3 

p'^Tj^_2P(x" + hf“) 

T -n 

= P T, ,PU k = 

k-3 

— n 

-X 

k = 

Tvn 


^ -OT^U 

k = 


Description of the Algorithm 

Calculation of the partial derivatives is performed simultaneously with the 
LSOM coordinate generation in the following order. 


Compute matrices C,B, and D, given by (15), (16), and (20) 
Compute matrices 

R 


J = 




[ujl (1-W^ ) 


J = JP^ 
o 


3. 

4. 


-n 


Compute vector Z = Pu 
Form five vectors 

- - (0, ’'g2®2“^s3®3’ ^s2®3'’'’'s3®2^ 

= rh (8^,81X32.85^33) 

r 

T _ -u / n n„ V 

J3 h (g2r,gj^X^2> §1^2^ 


= (0,-Z2»Z2) 
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= (Z^,0,-Z^) 

= (U2,-U^,0) 

Here r = || . 

5. Compute 

(k =1,2,. ..,5) 

w = -j3P 
6 

y y ” JJ y 

Here ]x^ is the 2x1 matrix of the partial derivatives with respect to 5 ^ and thus, 

11= (y^,y2,ii3>y^,y5,yg,yy)- 


Note, that SCD/GCD calculations require only the first six pair of the PD. The 
partial derivatives with respect to time,y^, will be used only to generate the 
backward scan grid for Thematic Mapper. 

The nominal SCD and PD are computed in double precision and stored in single 
precision. Because the PD are changing very slowly over a frame, they may be 
computed -on a sparse grid followed by linear interpolation onto a finer grid. 

For instance, implementation of our technique for MSS requires calculation of PD 
on a 3x5 grid. 


THE NOMINAL COORDINATES AND THE PARTIAL DERIVATIVES 
FOR BACKWARD SCANS OF TM 

It should be remembered that application of the developed technique to 

Thematic Mapper data requires two sets of the nominal SCD and PD - for forward and 

backward scans. But actually only one set must be obtained by the direct lookpoint 

calculation: LSOM coordinates for, say, forward scans may be easily converted to 

LSOM coordinates for backward scans. Our calculations show also that, for sensor’s 

o 

misalignments less than .1 , the derivatives are practically same for both grids; 
for bigger misalignment the second set of the derivatives can be obtained by the 
linear interpolation of the first one. 


Let X (t.,) and X (t^) be LSOM coordinates for adjacent forward and backward 


scans at 


t?. 


me^t 


and ^^y^ves^pect±^7ely 


Note, that for the TM 


dt = - 


^ 2 


T ^ + T , = .132205 sec 
act round 


So, we will neglect changes in the attitude angles during fit. 
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Considering Ax^ - X (t 2 )-X_(tj^) as a function of changes in the spacecraft position, 
sensor pointing, ana effects of the earth rotation, we may represent it as 




’4 "1 


1 1 


2 1 

+ •‘i®p + >'2 « 


A A 

Here 9 and 9 are fictitious pitch and roll angles, reflecting a difference in 
sensor s pointing at tj^ and t 2 , and is the average orbital rate during 6t. 


Here we will denote the nominal pointing vector g” at moments of time ti and to 


as P_and q, respectively 
plane, ( 0 ,P 2 ,P 3 > 

cos 0 = ~ 

p 


The angle between their projections onto the 
and (0,q^,q ), can be written as 

r> ^ I -n ^ 


^ 2^2 ^ 3*^3 


(P^ f P^)^ 


(q^ 


+ 93 ) 


or, choosing the appropriate sign. 


0 - sin 0 = 

P P 


^3^2 ^2^3 

9 1^ 9 ^ 

(l-Pp^(l-qp' 


Analogously, 0^ may be expressed as the angle between projections of P and q onto 
the plane: 


A 

0 


- Sin 0 




For zero sensor’s misalignments 

2 P P 
2 3 


P 


(i-pj)^ 


A 

0=0 

r 


CONVERSION TO BASIC MAP PROJECTIONS 


It should be remembered, that completely corrected imagery eventually must be 
presented in two basic map projections, SOM and either UTM or PS, To generate 
correction data in a basic map projection, it is required to invert LSGM coordinates 
to geodetic latitude and longitude and then perform the standard mapping into 
desirable projection. 
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Noting, that the normalized look point vector, V, can be expressed through the 
geocentric latitude, (j>, and the longitude, \ as 

= cosX cosij) 

^2 ~ sinX cos<j) 

= s±u<\) , 

A. 

and, employing well known formula for the geodetic latitude ^ 

^ 2 -2 
tan $ “ tan 


one can obtain 


X = arctan 

(j> = arctan |a^ b ^ V^(l-V^) 
I e e 3 3 




For a given X^, V is computed by the inverted formulae (10) and (9). 


NUMERICAL RESULTS 
The Accuracy of the Method 

To evaluate the methods accuracy, differences between LSOM coordinates, computed 
directly on points of a grid, and those, corrected according to Eq.(l), were 
calculated for various spacecraft position and attitude deviations. It is convenient 
to characterize the upper level of errors by the maximal along-track (AT) and 
cross-track (CT) errors, which coorespond to the errors at the worst points of a 
frame (possible different for AT and CT errors) . It should be noticed, that the 
maximal errors always appear near the corner points and similar for TM and MSS. 

They are linearly dependent upon magnitude of deviations and practically independent 
upon WRS latitude. 

The actual position and attitude departures for Landsat-D are expected to be 
Ol^(a) for the pitch, roll, and yaw and less than 5km in the along and cross track 
directions. The radial departure is determined chiefly by the orbit fluctuations 
and it will not exceed 9.5km. Modeling shows, that for 6 j, - 62 “ ^3 “ 5 

5 ^ 9.5km, the corresponding maximal CT and AT errors have the 
order of 5m (CT = 4.87m, AT = 5.03m for MSS and CT = 4.97m, AT = 5.07m for TM) . 

The inversion from LSOM to geodetic coordinates produces insignificant additional 
errors, therby preserving the same order of errors in UTM and PS. projections. 

For TM, the forward to backward scan conversion results in CT and AT errors , 
less than .03m for zero sensor’s misalignments; for extremely large misalignments 
of the order of .5 , the maximal CT errors increase up to .5m. 
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Currently the SCD/GCD generation accuracy for Landsat-D are defined in terms 
of the^ average mean-squared errors (Im for TM and 1,5m for MSS), Table 1 represents 
the 9Q/o maximal errors and the standard deviations of the average errors for Thematic 
Mapper, obtained by stochastic modeling. Here the attitude angles were normally 
distributed with zero means and a = .01 , Two cases of radial deviations were 
considered: a constant equal to 9.5 km, and a more plausible value from a uniform 

distribution (-9.5, 9,5) km. Because the errors do not depend significantly upon 
distribution of the cross and along track deviations, the latter were kept constant 
at 5 km, 400 samples were used to establish results for each case. The table 
also represents a case when PD were computed on a 3x7 grid and then interpolated to 
a finer grid. The nominal SCD and PD for backward scans were recomputed from the 
data for forward scans. Note, that in all cases the standard deviations of the 
average errors are less 1 m and thus, the nominal SCD and PD, precomputed for mean 
orbit radius at the frame center, provide the geometric correction with the required 
accuracy. 


Timing 


On the VAX, the direct lookpoint calculations take about 11 msec per grid 
point, interpolation of PD - 1 msec, the nominal SCD to SCD/GCD correction- less 
than .5 msec, and Inversion from LSOM to geodetic coordinates - 1.1 msec per point. 
Application of our technique for MSS requires interpolation PD to a finer grid, 
two corrections in LSOM coordinates, and inversion to geodetic coordinates; 
altogether it takes about 3,1 msec per point. The direct on-line SCD and GCD 
calculation takes about 22 msec per point. 

It should be noted, that mapping to the SOM requires about 15 msec per point, 
which is considered excessive for on-line processing. This time may be significantly 
reduced if we take into consideration the fact that the LSOM closely approximates 
true SOM distances between points within each frame. The errors of the approximation 
are relatively small (less than 5m) and sufficiently regular to permit linear 
interpolation LSOM to SOM coordinates. It may be done by using a 9x9 grid of 
corrections, precomputed and stored in the Data Base (Ref. 1). 


CONCLUSIONS 

The SCD/GCD calculation technique is based on presentation of image 
distortions as a sum of nominal distortions and linear effects, caused by variation 
of the spacecraft position and attitude variables from their nominals. The 
implementation requires generation and storage of the nominal SCD and twelve (for 
MSS), or fourteen (for TM) matrices of PD for each distinct latitude of WRS, along 
one path. The maximal errors of the method do not exceed 5.1m at the worst point 
of a frame. The standard deviations of the average errors are less than Im. 

The speed of the processing and the accuracy that is achieved by this technique 
makes it an elegant solution in the production environment. 
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Table 1. 


The 90% maximal errors and the standard deviations of the average errors for 
constant and uniformally distributed radial deviations. 


Distribution 
of the radial 
deviation 

Interpolation 
of PD 

Forward sea 

ns 

Backward scans 

90% max 
errors (m) 

STD 

(m) 

90% max 
errors (m) 

STD 

(m) 

CT 

AT 

CT 

AT 

CT 

AT 



CT 

AT 

constant 

no 

2.76 



.49 

2.78 


.61 

.49 

yes 

3.09 

2.69 

.63 

.67 

3.11 

2.70 

.70 

.72 

uniform 

no 

1.31 

1.78 

.29 

.32 

1.33 

1.78 

.58 

.58 

yes 

2.19 

1.41 

.42 

.35 

2.19 

1.40 

.53 

.48 
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